In [1]:
import plotly as plt
plt.offline.init_notebook_mode()
#!jupyter nbconvert --to html --execute 'drive/My Drive/Colab Notebooks/Exercicios_Capitulo_02.ipynb' --ExecutePreprocessor.iopub_timeout=3000
In [2]:
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
FILE_PATH = 'datasets/ipec.csv'
#FILE_PATH = 'drive/My Drive/Colab Notebooks/datasets/ipec.csv'
#!pip install lifelines
In [3]:
# Fuction to generate survival time of a group of individuals
from lifelines.plotting import plot_lifetimes
def survival_time_plot(p_entry_individual, p_actual_lifetimes, p_current_time, p_censured, p_time_unit, p_person_type, p_title, p_width, p_height):
observed_lifetimes = np.minimum(p_actual_lifetimes, p_current_time)
death_observed = p_actual_lifetimes < p_current_time
fig = plt.figure(figsize=(p_width,p_height))
ax = plot_lifetimes(durations = observed_lifetimes,
event_observed = death_observed,
entry = p_entry_individual,
left_truncated = True,
sort_by_duration = False,
event_observed_color = "black",
event_censored_color = "r",
)
total_individuals = len(p_actual_lifetimes)
total_period = max(p_actual_lifetimes+p_entry_individual)
ax.vlines(p_current_time, 0, total_individuals, lw=1, linestyles=':')
ax.set_xlim(-0.5, total_period + 1)
ax.set_xticks(np.array(range(0, total_period + 1)))
ax.set_yticks(np.array(range(0, total_individuals)))
ax.set_xlabel(p_time_unit, fontsize=12)
ax.set_ylabel(p_person_type, fontsize=12)
ax.set_title(p_title, fontsize=16)
ax.set_yticklabels(np.array(range(0, total_individuals))+1)
for iterator in range(0, total_individuals):
entry = p_entry_individual[iterator]
lifetime = observed_lifetimes[iterator]
period = entry + lifetime
ax.text(entry + ((lifetime)/2).astype(int), iterator + 0.2, "T=%d" % lifetime, fontsize=10)
marker = 'x' if p_censured[iterator] == 1 else 'o'
ax.plot(period, iterator, fillstyle='none', **dict(color='black', marker=marker, markersize=10))
ax.vlines(total_period, 0, total_individuals, lw=1, linestyles='--')
ax.plot(0,0, label='X - Evento', color='white')
ax.plot(0,0, label='O - Censura', color='white')
ax.legend(frameon=False)
ax.grid(False)
print("Tempo de vida observado: \n", observed_lifetimes, "\n")
In [4]:
# Função que gera gráfico de tragetória com variáveis N(t) e Y(t) do processo de contagem
def individual_trajectory_plot(p_entry_individual, p_lifetimes, p_censured, p_time_unit):
total_individuals = len(p_lifetimes)
total_period = max(p_lifetimes + p_entry_individual)
for individual in range(0,total_individuals):
entrytime = p_entry_individual[individual]
lifetime = p_lifetimes[individual]
censured = 'none' if p_censured[individual] == 0 else 'black'
point = entrytime + lifetime
print("Indivíduo: %d, Trajetória de Sobrevivência: %d." %(individual+1, lifetime))
fig, ax = plt.subplots(figsize=(6,3))
# Y(t)
ax.plot([0, entrytime],[.5, .5], color='r', linewidth=1)
ax.plot([entrytime, point],[1, 1], color='b', linewidth=1)
if point < total_period:
ax.plot([point, total_period],[.5, .5], color='r', linewidth=1)
ax.scatter(point, 1, marker="o", color='black', facecolors = censured)
if entrytime > 0:
ax.vlines(entrytime, 0.5, 1, lw=1, linestyles='--')
ax.scatter(entrytime, 0.5, marker="o", color='black')
if point < total_period:
ax.vlines(point, 0.5, 1, lw=1, linestyles='--')
# N(t)
if censured == 'black':
ax.plot([0, point],[-1, -1], color='r', linewidth=1)
ax.plot([point, total_period],[-0.5, -0.5], color='b', linewidth=1)
if point < total_period:
ax.vlines(point, -1, -0.5, lw=1, linestyles='--')
ax.scatter(point, -0.5, marker="o", color='black')
if point < total_period:
ax.text(point+0.3, -0.7, "dN(t)", fontsize=14)
else:
ax.text((total_period/2), -0.7, "dN(t)=0", fontsize=14)
else:
ax.plot([0, total_period],[-1, -1], color='r', linewidth=1)
ax.text((total_period/2), -0.7, "dN(t)=0", fontsize=14)
# Zero line
ax.plot([0, total_period],[0, 0], color='black')
ax.set_yticklabels(["", "0", "1","","0","1","",""],rotation=90)
ax.set_xticklabels([])
ax.text(0, -1.3, 0, fontsize=10)
ax.text(entrytime, -1.3, entrytime, fontsize=10)
ax.text(point, -1.3, point, fontsize=10)
ax.text(total_period, -1.3, total_period, fontsize=10)
ax.set_xlabel(p_time_unit, fontsize=14)
ax.set_xlim(-0.5, total_period + 0.5)
ax.xaxis.set_label_coords(0.5,-0.15)
ax.set_ylabel(" $N_{%d}(t)$ $Y_{%d}(t)$ " %(individual+1,individual+1), fontsize=14)
ax.set_title("$I_{%d}$ -> $T_{%d}=%d$" %(individual+1,individual+1,point), fontsize=14)
plt.show()
Exercício 2.1 ¶
O tempo de aleitamento, isto é, o tempo decorrido desde o nascimento até o desmame, pode ser considerado como uma variável tempo de sobrevivência. Suponha que o tempo até o desmame, em meses, tenha sido registrado para 15 crianças e que todas entraram no estudo no mesmo dia:
Tempos de aleitamento: 6 12 10 3 5 1 6 8 1 5 2 2 5 8 1
Considerando que não houve censura:
- a) Represente graficamente os tempos de observação das 15 crianças.
- b) Represente as trajetórias dos primeiros cinco indivíduos utilizando as variáveis N(t) e Y (t) do processo de contagem.
- c) Como você construiria um banco de dados para analisar estes dados pelo processo clássico?
- d)E pelo processo de contagem?
a) Represente graficamente os tempos de observação das 15 crianças.¶
In [5]:
# Dados de sobrevivência
entry_individual = np.zeros(15, dtype=np.int)
actual_lifetimes = np.array([6,12,10,3,5,1,6,8,1,5,2,2,5,8,1])
censured = np.ones(15, dtype=np.int)
current_time = max(actual_lifetimes) + 1 #This variable defines the time window
# Plot survival time
survival_time_plot(entry_individual, actual_lifetimes, current_time, censured,
"Meses", "Crianças", "Tempo de Aleitamento",
8, 6) # Plot size
b) Represente as trajetórias dos primeiros cinco indivíduos utilizando as variáveis N(t) e Y (t) do processo de contagem.¶
In [6]:
# Plot individual trajectory
individual_trajectory_plot(entry_individual, actual_lifetimes, censured, "Meses")
c) Como você construiria um banco de dados para analisar estes dados pelo processo clássico?¶
In [7]:
df = pd.DataFrame(np.array([6,12,10,3,5,1,6,8,1,5,2,2,5,8,1]))
df = df.reset_index()
df = df.rename(columns={0: 'Tempo', 'index':'Criança'})
df['Status'] = 1
df['Criança'] = df['Criança'] + 1
df
Out[7]:
d)E pelo processo de contagem?¶
In [8]:
df = pd.DataFrame(np.array([6,12,10,3,5,1,6,8,1,5,2,2,5,8,1]))
df['Tempo Inicial'] = 0
df = df[['Tempo Inicial', 0]]
df = df.reset_index()
df = df.rename(columns={0: 'Tempo Final', 'index':'Criança'})
df['Status'] = 1
df['Criança'] = df['Criança'] + 1
df
Out[8]:
Exercício 2.2 ¶
Considere agora o tempo de sobrevivência de 15 pacientes submetidos à hemodiálise. Neste caso, a variável de interesse é o tempo desde a primeira diálise até o óbito (em meses).
2, 4, 29+, 6+, 3, 1, 1, 2, 3, 9+, 10, 11, 5+, 5, 1
Os tempos censurados (censura a direita) estão indicados pelo sinal de +.
Considere que todos os pacientes entraram juntos no início do estudo.
- a) Represente graficamente os tempos de observação dos pacientes, utilizando a notação clássica.
- b) Represente as trajetórias dos primeiros cinco indivíduos utilizando as variáveis N(t) e Y (t) do processo de contagem.
- c) Construa um banco de dados para analisar esses dados pelo processo clássico.
- d) Reconstrua o banco de dados, na formas de processo de contagem, considerando que esta é uma coorte aberta e as datas de início de observação de cada paciente são dadas na tabela abaixo
a) Represente graficamente os tempos de observação dos pacientes, utilizando a notação clássica.¶
In [9]:
# Dados de sobrevivência
entry_individual = np.zeros(15, dtype=np.int)
actual_lifetimes = np.array([2, 4, 29, 6, 3, 1, 1, 2, 3, 9, 10, 11, 5, 5, 1])
censured = np.array([1,1,0,0,1,1,1,1,1,0,1,1,0,1,1])
current_time = max(actual_lifetimes + entry_individual) + 1 #This variable defines the time window
# Plot survival time
survival_time_plot(entry_individual, actual_lifetimes, current_time, censured,
"Meses", "Pacientes", "Óbito",
10, 6) # Plot size
b) Represente as trajetórias dos primeiros cinco indivíduos utilizando as variáveis N(t) e Y (t) do processo de contagem.¶
In [10]:
# Plot individual trajectory
individual_trajectory_plot(entry_individual, actual_lifetimes, censured, "Meses")
c) Construa um banco de dados para analisar esses dados pelo processo clássico.¶
In [11]:
df = pd.DataFrame(np.array([2, 4, 29, 6, 3, 1, 1, 2, 3, 9, 10, 11, 5, 5, 1]))
df = df.reset_index()
df = df.rename(columns={0: 'Tempo', 'index':'Paciente'})
df['Status'] = np.array([1,1,0,0,1,1,1,1,1,0,1,1,0,1,1])
df['Paciente'] = df['Paciente'] + 1
df
Out[11]:
d) Reconstrua o banco de dados, na formas de processo de contagem, considerando que esta é uma coorte aberta e as datas de início de observação de cada paciente são dadas na tabela abaixo.¶
In [12]:
df = pd.DataFrame(np.array([2, 4, 29, 6, 3, 1, 1, 2, 3, 9, 10, 11, 5, 5, 1]))
df['Tempo Inicial'] = np.array([0, 1, 12, 3, 1, 7, 0, 1, 1, 3, 12, 4, 4, 1, 22])
df = df.reset_index()
df = df.rename(columns={0: 'Tempo', 'index':'Paciente'})
df['Paciente'] = df['Paciente'] + 1
df['Tempo Final'] = df['Tempo Inicial'] + df['Tempo']
df['Status'] = np.array([1,1,0,0,1,1,1,1,1,0,1,1,0,1,1])
df
Out[12]:
Exercício 2.3 ¶
Suponha que, em um hospital, 6 pacientes HIV positivo são acompanhados ao longo de um ano.
No gráfico abaixo, as linhas horizontais representam o tempo de acompanhamento de cada paciente. Represente, utilizando as variáveis N(t) e Y (t), as trajetórias de cada um dos 6 pacientes.
In [13]:
# Dados de sobrevivência
entry_individual = np.array([3,0,5,0,6,3])
actual_lifetimes = np.array([9,6,3,12,6,4])
censured = np.array([0,1,0,0,0,1])
current_time = max(actual_lifetimes) + 1 #This variable defines the time window
# Plot survival time
survival_time_plot(entry_individual, actual_lifetimes, current_time, censured,
"Meses", "Pacientes", "Óbito por HIV",
10, 6) # Plot size
In [14]:
# Plot individual trajectory
individual_trajectory_plot(entry_individual, actual_lifetimes, censured, "Meses")
In [15]:
entry_individual = np.array([2,0])
actual_lifetimes = np.array([2,12])
censured = np.array([1,0])
individual_trajectory_plot(entry_individual, actual_lifetimes, censured, "Meses")
print("Qual foi o mês de entrada do paciente no estudo?", 2)
print("Em que mês ocorreu o desfecho?", 4)
print("Quais os meses em que o paciente estava sob risco de óbito?", 2, 3, 4)
In [16]:
entry_individual = np.array([0,0])
actual_lifetimes = np.array([8,12])
censured = np.array([1,0])
individual_trajectory_plot(entry_individual, actual_lifetimes, censured, "Meses")
print("Qual foi o mês de entrada do paciente no estudo?", 0)
print("Em que mês ocorreu o desfecho?", 8)
Exercício 2.6 ¶
Considere os dados de 6 indivíduos em uma coorte para a avaliação de fatores de risco para a infecção por HIV.
O estudo iniciou-se em 01/06/2008 quando 2 indivíduos entraram no estudo, com resultados negativos para a infecção.
Um deles ao realizar o segundo exame, um mês depois, teve resultado positivo, o outro positivou depois de 2 meses.
Os outros 4 indivíduos realizaram seus primeiros exames em 01/07, 01/09, 01/08 e 01/08, e as suas respectivas datas de soroconversões foram 2, 3, 2 e 5 meses, após a entrada.
In [17]:
entry_individual = np.array([0,0,1,3,2,2])
actual_lifetimes = np.array([1,2,2,3,2,5])
censured = np.array([1,1,1,1,1,1])
current_time = max(actual_lifetimes+entry_individual) + 1 #This variable defines the time window
a) Construa um banco de dados no formato clássico.¶
In [18]:
df = pd.DataFrame(actual_lifetimes)
df = df.reset_index()
df = df.rename(columns={0: 'Tempo', 'index':'Indivíduo'})
df['Status'] = censured
df['Indivíduo'] = df['Indivíduo'] + 1
df
Out[18]:
b) Construa um banco de dados em formato de processo de contagem..¶
In [19]:
df = pd.DataFrame(actual_lifetimes)
df['Início'] = entry_individual
df = df[['Início', 0]]
df = df.reset_index()
df = df.rename(columns={0: 'Fim', 'index':'Indivíduo'})
df['Status'] = censured
df['Indivíduo'] = df['Indivíduo'] + 1
df['Fim'] = df['Fim'] + df['Início']
df
Out[19]:
c) Represente graficamente os tempos observados desses indivíduos.¶
In [20]:
# Plot survival time
survival_time_plot(entry_individual, actual_lifetimes, current_time, censured,
"Meses", "Pacientes", "Soropositivo HIV",
10, 6) # Plot size
d) Represente as trajetórias dos seis indivíduos utilizando as variáveis N(t) e Y (t).¶
In [21]:
individual_trajectory_plot(entry_individual, actual_lifetimes, censured, "Meses")
e) Quais são os indivíduos em risco no 5º mês do estudo?
4 e 6
f) Que tipo de censura ocorreu neste estudo?
Censura intervalar
g) Em que situação haveria truncamento?
Incluíndo indivíduos já HIV positivos, ainda vivos no momento do primeiro exame.
Load dataset¶
In [22]:
df = pd.read_csv(FILE_PATH, sep=';')
df.describe()
Out[22]:
In [23]:
df.head()
Out[23]:
Function: surv¶
In [24]:
# Imprime variável de responsta e devolve em um array
def surv(p_df, p_ini, p_fim, p_tempo, p_status):
tempDf = pd.DataFrame()
classic_list = []
count_list = []
# Classic format
if p_tempo != '':
result = ""
for index, row in p_df.iterrows():
value = ("%s" % row[p_tempo]) + ("+" if row[p_status] == 0 else "")
result = result + "\t" + value
classic_list.append(value)
tempDf['Classic'] = np.array(classic_list)
print(result)
# Count process format
line_counts = 0
if p_ini != '':
result = ""
for index, row in p_df.iterrows():
value = ("(%s,%s" % (row[p_ini],row[p_fim])) + ("+" if row[p_status] == 0 else "") + "]"
result = result + "\t" + value
count_list.append(value)
tempDf['Score'] = np.array(count_list)
print(result)
return tempDf
a) No formato clássico¶
In [25]:
dfAux = surv(df, '', 'fim', 'tempo', 'status')
b) No formato de processo de contagem¶
In [26]:
dfAux = surv(df, 'ini', 'fim', '', 'status')
c) Observe a saída nos dois formatos. Quais os tempos e status do 1o e 79o pacientes em cada um dos formatos?¶
In [27]:
dfAux = surv(df, 'ini', 'fim', 'tempo', 'status')
print("\nPaciente 1")
print("Tempo processo clássico: " + dfAux.Classic[0])
print("Tempo processo de contagem: " + dfAux.Score[0])
print("")
print("Paciente 79")
print("Tempo processo clássico: " + dfAux.Classic[78])
print("Tempo processo de contagem: " + dfAux.Score[78])
In [0]:
In [28]:
df = pd.read_csv(FILE_PATH, sep=';')
a) Verifique as dimensões dos dados¶
In [29]:
df.shape
Out[29]:
b) Veja quais são as variáveis que constam no arquivo de dados¶
In [30]:
df.columns
Out[30]:
In [31]:
df.head()
Out[31]:
In [32]:
for col in df.columns:
print(col)
print(df[col].unique())
c) Substitua as informações ignoradas codificadas com 9/99 ou I por NA¶
In [33]:
dfAux = df.copy()
dfAux.anotrat.replace({9:'NA'}, inplace=True)
print('anotrat: ')
print(dfAux.anotrat.unique())
dfAux.obito.replace({'I':'NA'}, inplace=True)
print('obito: ')
print(dfAux.obito.unique())
dfAux.escola.fillna(9, inplace=True)
print('escola: ')
print(dfAux.escola.unique())
dfAux.fillna('NA', inplace=True)
dfAux[:5]
Out[33]:
In [34]:
for index, row in dfAux[:10].iterrows():
print("O paciente %d é do sexo %s, tem %d anos e foi acompanhado por %d dias %s" %
(row.id, ("masculino" if row.sexo == 'M' else "feminino"),row.idade, row.tempo, ("até a data do seu óbito." if row.obito == 'S' else ".")))
d) Indique quais são as variáveis que estão erroneamente classificadas como numéricas (isso acontece quando números são usados para codificar as categorias).¶
In [35]:
dfAux = df.copy()
ipec_dict = {'status': {0:'censura',
1:'óbito'},
'escola': {0:'sem escolaridade', 1:'ensino fundamental',
2:'ensino médio', 3:'ensino superior'},
'risco': {0:'homossexual masculino', 1:'usuário de drogas injetáveis', 2:'transfusão',
3:'contato sexual com HIV+', 5:'hétero c/múltiplos parceiros', 6:'dois fatores de risco'},
'acompan': {0:'ambulatorial/hospital-dia', 1:'internação posterior', 2:'internação imediata'},
'tratam': {0:'nenhum', 1:'mono', 2:'combinada', 3:'potente'},
'doenca': {1:'pcp', 2:'pcp pulmonar', 3:'pcp disseminada', 4:'toxoplasmose', 5:'sarcoma',
7:'outra doença', 8:'candidíase', 9:'duas doenças', 10:'herpes', 99:'definido por cd4'},
'propcp': {0:'sem profilaxia', 2:'primária', 3:'secundária', 4:'ambas'},
'obito': {'S':'óbito', 'N':'não óbito', 'I':'ignorado'}
}
dfAux.anotrat.replace({9:'NA'}, inplace=True)
dfAux.fillna('NA', inplace=True)
dfAux.replace(ipec_dict, inplace=True)
dfAux.head()
Out[35]:
e) Calcule as medidas resumo das variáveis presentes no banco¶
In [36]:
from pandas.api.types import is_string_dtype
from pandas.api.types import is_numeric_dtype
for col in df.columns:
if is_string_dtype(dfAux[col]):
print('-----------------------------')
print('Feature: %s' % col)
print(dfAux[col].describe())
print('-- Counts')
print(dfAux[col].value_counts())
else:
print('-----------------------------')
print('Feature: %s' % col)
print(dfAux[col].describe())
f) Construa também tabelas de freqüência das variáveis de maior interesse¶
In [37]:
pd.DataFrame(dfAux.groupby('escola').size())
Out[37]:
In [38]:
pd.DataFrame(dfAux.groupby('tratam').size())
Out[38]:
In [39]:
pd.DataFrame(dfAux.groupby(['status']).size())
Out[39]:
In [40]:
pd.DataFrame(dfAux.groupby(['sexo','escola']).size())
Out[40]:
In [41]:
pd.DataFrame(dfAux.groupby(['risco','sexo']).size())
Out[41]:
g) Construa gráficos para o tempo de sobrevivência e para outras variáveis de interesse¶
In [42]:
import plotly.express as px
fig = px.box(df, x = 'sexo', y='idade', points='all')
fig.update_layout(title_text='Idade por gênero',
xaxis_title_text='Sexo',
yaxis_title_text='Idade'
)
fig.show()
In [43]:
fig = px.box(dfAux, x='sexo', y='idade', color='risco')
fig.update_layout(title_text='Idade por sexo e risco',
xaxis_title_text='Sexo',
yaxis_title_text='Idade'
)
fig.show()
In [44]:
fig = px.histogram(dfAux, x="obito", y="id", histfunc='count', facet_col='sexo', opacity=0.5)
fig.update_layout(title_text='Sobrevivência',
xaxis_title_text='Ocorrência', # xaxis label
yaxis_title_text='Quantidade', # yaxis label
bargap=0.05, # gap between bars of adjacent location coordinates
)
fig.show()
In [45]:
fig = px.histogram(dfAux, x="tempo", nbins=10, labels={'tempo':'Dias'}, opacity=0.5)
fig.update_layout(title_text='Sobrevivência',
xaxis_title_text='Dias', # xaxis label
yaxis_title_text='Frequência', # yaxis label
)
fig.show()
In [46]:
fig = px.histogram(dfAux, x="tempo", nbins=10, color='sexo', labels={'tempo':'Dias'}, opacity=0.75)
fig.update_layout(title_text='Sobrevivência',
xaxis_title_text='Dias', # xaxis label
yaxis_title_text='Frequência', # yaxis label
bargap=0.05, # gap between bars of adjacent location coordinates
)
fig.show()
In [47]:
fig = px.histogram(dfAux.loc[dfAux['status'] == 'óbito'], x="tempo", nbins=10, color='sexo', labels={'tempo':'Dias'}, opacity=0.75)
fig.update_layout(title_text='Eventos',
xaxis_title_text='Dias', # xaxis label
yaxis_title_text='Frequência', # yaxis label
bargap=0.05, # gap between bars of adjacent location coordinates
)
fig.show()
In [48]:
fig = px.histogram(dfAux.loc[dfAux['status'] == 'censura'], x="tempo", nbins=10, color='sexo', labels={'tempo':'Dias'}, opacity=0.75)
fig.update_layout(title_text='Censurados',
xaxis_title_text='Dias', # xaxis label
yaxis_title_text='Frequência', # yaxis label
bargap=0.05, # gap between bars of adjacent location coordinates
)
fig.show()